Fractal valence bond loops in a long-range Heisenberg model at criticality 
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We present a valence bond theory of the spin-S quantum Heisenberg model. For nonfrustracting, 
local exchange and dimension d > 1, it predicts a resonating ground state with bond amplitudes 
h{r) ~ (a^ + r^^^'Pi'^ and decay exponent p = d + 1. Different values of p can be achieved by 
introducing frustrating (p > d + 1) or nonfrustrating (p < d + 1) long-range interactions. For d — 2, 
but not d = 3, there is a critical value of the decay exponent pc above which the ground state is 
a spin liquid. The phase transition is analogous to quantum percolation, with fractal valence bond 
loops playing the role of percolating clusters. The critical exponents are continuously tunable along 
the phase boundary p = Pc(a, S). 



Singlet product states were introduced in the early 
days of quantum chemistry to describe valence bond- 
ing in molecules P, 0]. These "valence bond states" 
were also applied to translationally invariant systems and 
proved especially useful for studying interacting S = ^ 
quantum spin chains [1, 13| • A featureless resonating va- 
lence bond (RVB) wavefunction, consisting of a superpo- 
sition of all possible configurations of short range valence 
bonds, was later promoted by Anderson as a possible 
ground state for low-coordination-number antiferromag- 
nets [5, 6] and as apotential route to superconductivity 
in the cuprates 0, 

Liang, Doucot, and Anderson (LDA) extended the 
RVB picture to include bonds on all length scales by 
assuming that the weight associated with each config- 
uration could be factorized into a product of individual 
bond amplitudes^ach expressed as a single function of 
the bond length [9|. They considered various functional 
forms for the two-dimensional Heisenberg model and con- 
cluded that the optimal amplitudes fall off as a powerlaw: 
h(r) ^ r^P. Their variational calculation was not suffi- 
ciently accurate to resolve the correct exponent, which is 
almost certainly p = 3 



11|,|12|. 



In this Letter, we begin by providing a formal justifica- 
tion for the LDA ansatz. The calculation presented here 
is a mean field treatment of the Heisenberg model based 
on valence bond creation and annihilation operators [ist ; 
these obey an unusual operator algebra that reproduces 
the properties of the overcomplete valence bond basis. 
We show that a broad class of spin-S* models with two- 
spin interactions have ground state wavefuctions that are 
very close to being factorizable-amplitude RVB states. 
Antiferromagnetically ordered states attain perfect fac- 
torizability in the large spin limit. 

For models with local, nonfrustrating interactions, the 
mean field equations can be solved analytically. The am- 
plitudes are generically radially symmetric and of the 
form h{r) — (a^ -|- r'^)~P^^, where the core size of the 
distribution is a = 1/ ^/d and the powerlaw tail has ex- 
ponent p = d+ 1. Here d is the dimension of the lattice. 
The only exception is for d = 1, in which case either 



the bond lengths are gaussian distributed (even 25') or 
beyond mean field prediction (odd 2S). In summary. 
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d=l,2S even (1) 
d = 2 (2) 

d = 3 (3) 



The ^ appearing in Eq. ^ is the spin correlation length 



^2^^^ J_g.(2S+l)^ 



(4) 



We then ask the question, are there couplings that 
preserve the h{r) = {a'^ + r^)^P^^ form but accommodate 
arbitrary values of p7 Indeed, it turns out that tuning p 
away from d-l-l is equivalent to introducing a long-range, 
bipartite interaction j in opposite sublattices) 

^.=7. + (l-7.) ^^^^^^ + ;'"\ (5) 

where 7^ is the nearest neighbour matrix and a{p),f3(p) 
are smooth, positive functions. The beyond-nearest- 
neighbour part is frustrating when Sp = p^d~l>0. 
For the square lattice, a(3 + 5p) — 4 + O.SASp and 
(3{3 + 6p) =0.7- 0.15 (5p. 

The function h{r) — {a^ + r^)~P^^ characterizes a fam- 
ily of RVB states interpolating between the classical Neel 
state {p 0) and Anderson's short-range RVB state 
(jj — > 00). For d — 2, these two limits are separated 
by a phase transition: there is a critical line p = pda, S) 
beyond which the system is quantum disordered. See 
Fig. [TJ Along the critical line, the valence bond loops 
are scale invariant and have fractal dimension Df. Their 
length distribution is characterized by entropy and loop 
tension exponents, r(a, S) and CT(a, S), in terms of which 
all other critical exponents can be expressed. In higher 
dimensions, where even the short-range RVB state is Neel 
ordered, there is no transition. 
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FIG. 1: A phase diagram for the RVB wavefunction with 
ampUtudes h{r) — (a^ + r^)~P^'^ on the square lattice, r 
is the length of a valence bond connecting spins in opposite 
sublattices. The phase boundary is drawn for S — ^ (solid 
line) and for S = 1, |, 2 (dotted lines, from left to right). The 
location of the ground state (x) can be adjusted horizontally 
by perturbing [as per Eq. ([5])] about the nearest-neighbour 
Heisenberg model. Points at which the critical exponents have 
been measured (see Table are marked (•). 



Mean field — The valence bond creation and annihila- 
tion operators XiJ and xfj > discussed in Ref . [S, form an 
overcomplete basis for any system consisting of an even 
number of SU(2) spins. The operator XiJ — iXijiXlj) 
creates a singlet (/i — 0) or triplet — 1,2,3) pair at 
sites i and j out of the spinless vaccuum. More precisely, 



X^^lvac) 



(6) 



where r'^ — (icr^, icr^, 1, — icr^) is a vector of 2 x 2 matrices 
related to the Pauli matrices by t° — icr^ and r — cra^. 

Valence bond operators with no indices in common 
commute. Those with matching site indices obey a 
bosonic commutation relation, whereas those with only 
one index in common obey a more complicated rule: 



[Xij : Xij 1 



r fj, pt »i 



vac) 



vac) 



(5^" I vac 
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(7) 

(8) 



Here, T•'P^''^ ^ ^trT^^TPr^'^T 



takes the values 0,±1. 
Equation ([5]) encodes the overcompleteness of the basis. 

Valence bond operators can be used to construct 
higher-spin states, so long as there are exactly 2S bonds 
emerging from each site: J^fi Xij Xij = 25*. Then each 
valence bond state is characterized by a set v specifying 
the endpoints and singlet/triplet character of each bond. 
That is, \v) — S[[\xij] |vac), where the product is over 
all {i,j,fj,) e V and S denotes a symmetrization of the 
operators with respect to all orderings 1^, 17 1. 



Since singlet bonds have a direction associated with 
them {x% = ~Xijj whereas Xji = Xij)i we restrict our 
discussion to bipartite lattices (A^ spins in each sublat- 
tice) and adopt the convention that all bonds originate 
in the A sublattice and terminate in the B sublattice. 
Hence, operator indices appear in the standard order XiJ 
with i G A and j € B, and the transform Xq is defined 
at only TV wavevectors in a reduced BriouUin zone. 

An arbitrary state can be expressed as a linear combi- 
nation (not unique) ji/)) = iIj{v)\v). Factorizability in 
the LDA sense means that the coefficient tp{v) = JJa'^j 
can be decomposed into a product of individual bond 
amplitudes a^-. Even under the assumption of perfect 
factorization, this wavefunction is nontrivial, since there 
are geometrical constraints implicit in tiling the lattice 
with hard-core dimers. 

Those constraints can be made explicit in the following 
way. The bond configuration sum can be achieved by 
exponentiation of a weighted one-bond operator: 
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XiVlvac). 



(9) 



The Gutzwiller projection operator G = HisA 6{Ni~2S) 
filters out all configurations that do not have exactly 25* 
bonds emerging from each lattice site. This is similar in 
spirit to the construction of Anderson's projected BCS 
wavefunction 0], but here the underlying operators are 
not paired fermions. The projection step can be imple- 
mented by introducing a gauge field at each bond end- 



point: X 
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and G = /i:»[0]e-'^^ ^^s^,^ 



Equation ([9]) can be written compactly as — G\a) 
with I a) = exp^q ^ a(^Xq^|vac). If we relax G and en- 
force the number constraint on average (fixing the overall 
bond number at 2SN), then \a) can be viewed as a mean 
field approximation to Since Ni ^ d/d4>i, the total 
bond number operator N — J^ieA related to the 

overall normalization of the bond amplitudes: 



N\a) = ^\ua) 



E 



'-X^^la). 



qA,q 



(10) 



Let us assume that the spin-spin interactions in the 
Hamiltonian act only between sites in opposite sublat- 
tices: H = J^jS^ ■ Sj - X{N - TV). We have in- 
troduced a Lagrange multiplier A coupled to the shift in 
bond number. In the state |a), the valence bond opera- 
tor has an anomalous expectation value = Jy(xf,). 



For a singlet ground state, the various triplet components 



vanish: Af^ 



JSSt"'" and 



The 



expectation value of Si ■ Sj between valence bond states 
depends on whether i and j are in the same valence bond 
loop 0. At the mean field level, 



C., = (-1)^+^(8,. S,) 



(5+1)%, (11) 
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where X„ = f[h^/{l - i|/iqP)] and X,k = f[V2/{l - 
i|/iqp)] (with i,k ^ A and j G B) measure the probabil- 
ity that two sites are connected by a chain of bonds. By 
the variational principle, d{H)/dh* — yields 



(X°) = - 
and d{H) /dX 



A 







c.q= JA2-2|Aop, (12) 



yields 
1 ^ A 



1 + 2S'. 



(13) 



If Jij ~ Jjij , the dispersion is = (e^ + 1 



>\l/2 
iq) ^ 

in units where |Aq^g| = l/\/2 and = 1 + e^. In 
d = 1, there is an important difference between integer 
and half-integer values of 5* For even 25", excita- 

tions are gapped and Eq. (fTS]) reads — ^ log ^ = 2S' -I- 1, 
which leads to Eq. ([T|). For odd 25', the mean field 
breaks down; the system is gapless [l9j, but LOq ^ q 
causes the bond constraint equation to diverge. [The 



bosonization result, 



1/2 



(logr; 



sug- 



gests h{r) ^ r^'^/^(loglogr-|~0(l))]. For d > 2, however, 
the gap e = OlL^"^) vanishes faster than the wavevec- 
tor spacing n/L. Thus, as L ^ oo, the divergence in 
the bond number equation is confined to the zero mode 
and the stucture of the equation mimics that of so-called 
"sublattice-symmetric spinwave theory" jl^ : 



q 



+ c=l + 2S. 



(14) 



The details of the lattice enter as a single geometrical 
constant, c = 1.393 (square), c = 1.1567 (cubic), etc. 
The staggered moment extracted from the large separa- 
tion limit of Xij = (xij) is M = S — ^(c — 1), the usual 
spinwave result [l^. e — *■ leads to the long- wavelength 
behaviour log(/iq/-\/2) = —q/Vd+0{q^), with deviations 
from radial symmetry not appearing until O(g^). The 
corresponding real-space bond amplitudes for c? = 2, 3 
[Eqs. ([2]) and ([3])] are found by computing the Hankel 
transform. In general, the presence of gapless, linearly 
dispersive excitations implies h{r) ~ l/r''+i. The core 
size a = 1/Vd = -\/2ws/| A^^qI is proportional to the 
spin- wave velocity Vg. 

We now consider the broader class of wavefunctions 
with amplitudes h{r) = (a^ -I- r^) and compute the 
corresponding p-a phase diagram. In the d ~ 2 case 
(Fig. [T]), there is a quantum disordered phase at suffi- 
ciently large p for all values of a and S. The critical ex- 
ponents (/?, V, z can be measured directly; see Fig. [2]) vary 
continuously and monotonically along the line of transi- 
tions with /3/0.8, 1^X1.25, and z\1.36 as a\0. The 
loss of long range order can be attributed to a qualitative 
change in the distribution of valence bond loops — a van- 
ishing number of system-spanning loops on the magnetic 
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FIG. 2: (Top-left) The staggered moment M and spin gap 
ajq=o on a square lattice of size L = 192 are plotted for 25 = 1, 
a — l/\/2. The critical value pc ~ 3.6256(5) is determined 
from crossing points of the ratio R = X{^, ^)/ X{j , j); 
the curves L = 64, 96, 128, 192, 256 are shown in the inset. 
(Bottom-left) Critical exponents v (□), /3 (O), and z (A) 
for 25 = 1 are plotted alongside the predicted values (solid 
lines) of the exponents based on a [l,l]-Pade approximant fit 
of r and a (dashed lines). (Right) The staggered moment 
on the L = 32 cubic lattice is shown for 2S — 1,2,3,4 and 
2a = 0, 1, . . . , 5. In the limit L,p, 1/a — > oo, its value bot- 
toms out at Moo = S — 0.52, indicated by dotted lines and an 
arrow. 



side becomes a macroscopic number of small loops on the 
spin liquid side [l3| . At criticality, the problem is for- 
mally equivalent to that of percolation, with self-similar 
valence bond loops of fractal dimension Dj = d—fi/v (20| 
playing the role of percolating clusters (although loops, 
unlike clusters, are all backbone, i.e., free of dangling 
spins [13]); fluctuations of the loop gas substitute for 
the geometric average over disorder configurations. More 
precisely, the analogy is with quantum, rather than clas- 
sical percolation: the bare dynamical exponent zq = 1 
should be renormalized according to z = 4>Df = Df, 
where 4> = zo/(2 — ^o) = 1 [l^- This is borne out by 
measurements of the critical exponents, listed in Table [H 

Near criticality, the loops are distributed according 
to 71; ~ Z^'^exp[— (pc — py^'^l], ni being the number 
of loops of length I. The distribution is parameterized 
by an entropy exponent t and a loop tension exponent 
cr, which are related to the conventional exponents by 
J/ = (r - l)/da, I3^[t- 2)/cr, and z = Df ^ d/{T - 1). 
The bottom- left panel of Fig. [5] illustrates the excellent 
agreement between measured exponents and their pre- 
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TABLE I: Critical exponents — S — ^, square lattice 



a 


Pc 


u 


/3 


z 


Df 





3.3677(2) 


1.25(4) 


0.80(2) 


1.361(7) 


1.36(3) 


I/V2 


3.6256(5) 


1.7(1) 


0.69(4) 


1.594(6) 


1.59(3) 


1 


3.800(1) 


2.0(1) 


0.56(3) 


1.718(6) 


1.72(2) 


2 


4.24(1) 


2.6(1) 


0.26(1) 


1.902(4) 


1.90(5) 




p = 3.7 -p — h p = oo 



dieted values based on a [l,l]-Pade fitting form for r and 
a. Note that the anomalous dimension, determined from 
the relation t] = 2f3/iy— z + 2 — d, is strictly negative and 
goes as 77 ss —a. 

On the disordered side of the transition, the system is a 
gapped spin liquid. In the p ^ 00 limit, the short-range 
RVB state on the square lattice has a spin correlation 
function 

^5(5 + 1)^^, (15) 

where ^(i) = 0.83, ^(1) = 2.0, C(|) = 4.4, and log^iS) - 
1.5 S for large spin. In contrast, long-range antiferromag- 
netic order survives in d = 3 for all p. This is true at the 
mean field level provided that S > ^ (see Fig.[2|); a direct 
evaluation of the short-range RVB ground state (using an 
efficient valence bond worm algorithm) confirms that the 
S = ^ case is also Neel ordered. 

We now run the mean field equations in reverse and 
solve, via 

A l + i|^P' A F[(x^)] ' ^ ' 

for the interaction that stabilizes ground states with 
a — and arbitrary p (parameterizing, e.g., the 

horizontal trajectory in Fig. [T]). The result is the in- 
teraction given in Eq. Numerical estimates on the 
square lattice show that the long-range part goes as 
[0.83(4 - a) + 0.21(4 - + ■ ■ ■ ]/rfj. The solutions are 
characterized by o^q ^ l/(Xq) ~ 9^ and 1 — ^ q^^ , 
where 9 = 1 -|- (a — 4) + 0{a — 4)^ « a — 3. For nonfrus- 
trating interactions {a < 4), the dispersion is sublinear 
and for frustrating interactions (a > 4), it is superlinear. 
As the excitations soften, there is a point {d — 2 only) at 
which the zero mode of the constraint equation ceases to 
hold macroscopic weight. Here, the magnetic order van- 
ishes and a gap opens in the excitation spectrum. For 
5 = i, that point occurs at 9c == 1.5 and at large S, 
9c-2~j^. Oddly, Eq. © (with J,, ^ ~3.6/rff) 
takes us up to but not through this critical point. Entry 
into the p > Pc phase corresponds to the development 
of an unusual concentric-ringed, cloverleaf pattern of al- 
ternating frustrating and nonfrustrating interactions, as 
shown in Fig. [3l Other Hamiltonians with tunable inter- 
action strengths potentially lead to different trajectories 
in the p~a plane. 



FIG. 3: The sign of the exchange interaction Jij is plotted 
for a — I/V2 and various values of p > Pc- Dark dots denote 
Jij > and light dots Jij < 0. 

A key point is that the long-range, frustrating terms 
in Eq. ^ push p to higher values (toward the disordered 
phase) while preserving the radial symmetry and posi- 
tivity of h{r). This is generally not true of short-range 
frustrating interactions, which tend to break the radial 
symmetry {h{r) — s- h{r)) and, at large couplin g st rength, 
the Marshall sign rule (3 r such that h{r) < 0) Since 
reduction of the continuous rotational symmetry to the 
discrete rotational symmetry of the lattice is a precursor 
to crystalline bond ordering, preservation of the radial 
symmetry is important if the liquid state is not to be 
preempted by a bond ordered one. 
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